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1.  Introduction 


The  restoration  of  degraded  images  has  many  fields  of  appli¬ 
cation,  including  space  imagery  and  biomedical  images  [1-2] . 

Several  approaches  have  been  suggested  in  the  literature  using 
recursive  Kalman  filter  type  algorithms  [3-5]  ,  inverse  filtering 
methods,  and  two-dimensional  discrete  partial  difference  equations 
[6-7] .  The  Kalman  filter  algorithms  insist  on  the  assumption  of 
an  underlying  causal  representation  for  the  images  displayed  as 
a  state  space  model.  The  various  Kalman  filter  methods  differ 
from  one  another  basically  regarding  the  assumptions  that  they 
make  about  the  autocorrelation  function.  Much  attention  has  been 
paid  to  exponential  separable  correlation  functions  in  the  litera¬ 
ture.  Also,  the  optimality  with  respect  to  MMSE  claimed  in  Kalman 
filter  type  algorithms  is  valid  only  when  the  parameters  of  the 
underlying  RF  models  are  exactly  known.  In  practice,  the  parameters 
of  the  model  are  estimated  from  the  images,  and  substituted  for 
the  unknown  parameters.  Optimality,  with  respect  to  some  criterion, 
can  be  claimed  only  if  estimation  of  the  parameters  and  minimization 
of  the  mean  square  error  are  tackled  simultaneously  so  as  to 
minimize  the  criterion  function.  Detailed  discussions  on  inverse 

filtering  approaches  may  be  found  in  [1-2]  .  _ 

Recently  [6-7]  ,  stochastic  representations  of  digital  images  ~  — 
by  finite  difference  approximations  of  partial  differential 
equations  (PDE)  have  been  used  to  develop  restoration  algorithms. 
Corresponding  to  the  PDE  classification  of  hyperbolic,  parabolic, 
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and  elliptic  systems,  three  different  algorithms,  viz.,  the 
causal,  semicausal  and  noncausal  algorithms,  have  been  developed. 

By  making  appropriate  assumptions  regarding  the  boundary  pixels, 
the  original  observations  are  decomposed  into  two  parts  such  that 
the  covariance  matrix  corresponding  to  one  set  is  diagonalized  by 
sine  transforms  and  the  other  part  corresponds  to  the  boundary 
response.  The  important  feature  is  that  the  eigenfunction  expan¬ 
sion  of  the  covariance  matrix  in  terms  of  sine  functions  enables 
fast  implementation  of  the  restoration  filters. 

We  consider  the  use  of  spatial  autoregressive  models  for  the 
image  restoration  problem.  Within  the  family  of  spatial  autore¬ 
gressive  models,  there  are  two  nonequivalent  classes  of  RF  models, 
the  so-called  simultaneous  autoregressive  (SAR)  [8-10]  models  and 
the  conditional  Markov  (CM)  [10-13]  models.  These  models  charac¬ 
terize  the  statistical  dependency  of  a  pixel  on  its  reighbors  by 
a  linear  weighted  sum  of  neighboring  pixels  and  additive  noise, 
independent  or  correlated  depending  upon  whether  the  models  are 
simultaneous  or  conditional.  The  structure  of  these  models  can  be 
explained  as  follows:  suppose  we  are  given  an  observation  set 
(y(s),  s  =  (i,j)fcfi},  «={s  =  (i , j )  l<i,j<M},  where  y(s)  is  the 
observation  at  location  s.  Then  the  expectation  of  y(s),  condi¬ 
tional  on  all  yts^),  s^fcfifS^s,  is  only  a  function  of  the  obser¬ 
vations  at  positions  belonging  to  the  neighbor  set  for  the  condi¬ 
tional  models,  whereas  for  the  simultaneous  models,  the  conditional 
expectation  is  a  function  of  observations  at  positions  belonging 


to  the  neighbor  set  and  some  extra  neighbors.  For  instance,  the 
conditional  expectation  of  a  four  neighbor  simultaneous  model  with 
dependence  on  east,  west,  north  and  south  neighbors  is  a  function 
of  eight  nearest  neighbors  and  second  nearest  neighbors  on  the  east, 
west,  north  and  south.  These  two  classes  of  models  are  non¬ 
equivalent,  in  that  given  a  SAR  model,  an  equivalent  CM  model  can 
always  be  found  (usually  with  more  parameters) ;  however,  the 
converse  is  not  true.  For  instance,  there  is  no  equivalent  SAR 
model  for  the  CM  model  with  dependence  on  the  east,  west,  north 
and  south  neighbors.  In  this  paper,  we  are  primarily  concerned 
with  the  use  of  SAR  models  for  the  digital  image  restoration 
problem.  Similar  results  for  CM  models  will  be  given  later.  For 
some  results  concerning  the  use  of  conditional  models  in  image 
restoration,  the  reader  is  referred  to  [5,14]. 

For  a  finite  image,  the  neighbors  in  all  directions  are  not 
defined  for  boundary  pixels.  Some  assumptions  have  to  be  made 
regrading  the  distribution  of  boundary  pixels.  In  this  report,  we 
assume  that  the  given  finite  images  are  represented  on  a  torus 
lattice.  This  representation  on  torus  lattices  leads  to  covariance 
matrices  that  have  block  circulant  structures.  Since  block  circu- 
lant  matrices  are  diagonalized  by  Fourier  vectors,  FFT  computations 
can  be  used  for  the  implementation  of  restoration  filters. 

Circulant  approximations  to  Toeplitz  covariance  structures  have 
been  considered  in  the  image  restoration  literature.  But  the 


block  circulant  structure  that  arises  in  our  case  is  due  to  the 
underlying  RF  model.  Consequently,  the  covariance  matrices  are 
exactly  block  circulant. 

Using  the  specific  representation  mentioned  above  for  the 
images,  the  restoration  problem  is  then  posed  as  one  of  minimiz¬ 
ing  the  mean  square  error  between  the  original  and  the  restored 
image.  This  assumes  that  the  SAR  model  is  completely  specified, 
i.e.,  the  parameters  characterizing  the  SAR  model  are  known.  In 
practice,  however,  these  parameters  are  unknown  and  are  estimated 
from  the  images.  For  SAR  models  that  include  dependency  in  all 
directions,  the  Jacobian  of  the  transformation  matrix  from  the 
noisy  variates  to  the  observations  is  not  unity,  leading  to  a  log 
likelihood  function  that  is  nonquadratic  in  the  parameters.  Direct 
minimization  of  the  log  likelihood  function  using  "off  the  shelf" 
computer  programs  might  lead  to  a  slow  convergence  rate.  Conse¬ 
quently,  an  iterative  scheme  suggested  elsewhere  [15-16]  is  used 
for  estimating  the  parameters. 

The  organization  of  the  paper  is  as  follows:  In  Section  2,  we 
pose  the  restoration  problem  using  the  SAR  models  and  derive  its 
optimal  filter,  given  the  parameters  of  the  SAR  model.  A  brief 
discussion  is  given  regarding  the  estimation  of  unknown  parameters 
in  SAR  models.  Performance  bounds  for  the  restoration  algorithms 
are  also  derived.  Several  examples  of  restorations  are  given  in 
Section  3.  Finally,  discussion  is  given  in  Section  4. 


2 .  Restoration  Schemes  Using  the  Simulatnneous  Random  Field  Models 
2 . 1  Derivation  of  the  Optimal  Filter 

Consider  an  MxM  zero  mean,  undegraded  image  represented  by  a 
two-dimensional  array  of  gray  levels,  {y (s) ,s=(i, j)  Q} ,fi={ s= (i , j ) , 
l<i,j_<M}  where  y(s)  is  the  gray  level  of  the  cell  s.  Since  for  a 
finite  image,  some  of  the  neighbors  for  the  boundary  points  are 
not  defined,  the  image  is  assumed  to  be  folded  into  a  torus  so 
that  (2.1)  is  satisfied: 

y[  (s)  +  (i-^ ,  j  ^)  ]=y[{  (s)  + (i^-M, j^-M)  }  (mod  M+l,mod  M+l)  ]  (2.1) 

Assume  that  the  original  image  y(s)  obeys  the  SAR  model 

y (s) =  E  0.  -y(s+(i, j) )+/pw(s) ,  s  fi  (2.2) 

(i,j)€N  1,3 

In  (2.2),  {w(s),s€fi)is  a  sequence  of  independent  and  indentically 
distributed  random  variables  of  zero  mean  and  unit  variance  and 
N  denotes  the  neighbor  set,  i.e.,  the  set  of  neighbors  on  which 
the  observation  at  y(s)  is  dependent.  Equation  (2.2)  is  charac¬ 
terized  by  a  set  of  parameters  {0.  .  ,  (i, j)fcN,p)  .  Typically,  N 

1 1  3 

might  include  the  members  { (-1 , 0) , (0 , 1) , (1 , 0) , (0 , -1) }  corresponding 

to  nearest  neighbors  in  north,  east,  south,  and  west  directions, 

or  N  might  be  N={ (-1 , 0) , ( 1 , 1) , (0 , 1) , (1 , 0) , (0 ,-l) } .  One  of  the 

characteristics  of  the  model  in  (2.2)  is  that,  when  the  parameters 

corresponding  to  symmetric  neighbor  sets  are  interchanged,  the 

resulting  model  has  second  order  properties  identical  to  the  ori- 

2 

ginal.  Denoting  y  and  w  as  M  xl  vectors  of  lexicographically 

ordered  arrays  { y  ( * ) >  and  (w(*))»  (2.2)  can  be  written  as 

B ( 0 ) y  =  /Pw  (2.3) 

2  2 

where  B(0)  is  an  M  xM  block  circulant  matrix  and  has  the 


following  structure: 


i 


where  each  of  the  component  matrices  is  circulant.  For  the  case 
N={  (0,1)  ,  (1,1)  ,  (0,-1)  , (-1,0)}  , 

1=circulant  (l,eQ  -1* 

B1  2=circulant  (9-j^  lf0,...,0) 

B^  M=circulant  (9_1  Q,0,...,0) 

and 

B,  -=0  j*l,2,...,M. 

f  J 

Assuming  that  B(0)  has  an  inverse,  we  have 

y  =  /p  B(9)_1w  (2.4) 


(2.6) 


Jit. 


T 

E (nn  )  =  C  (a  block  circulant  matrix) 
we  have 

x  =  HY+n  (2.7) 

We  assume  for  the  present  that  the  PSF  matrix  H  and  the  parameters 
(0,p)  in  (2.2)  as  well  as  the  matrix  C  in  (2.6)  are  completely 
known.  We  also  assume  that  the  degraded  image  fits  completely 
inside  their  recording  frame.  In  other  words,  the  problems  due 
to  truncation  of  the  degraded  image  by  a  finite  recording  frame 
are  not  considered.  The  restoration  problem  is  posed  as  follows: 
determine  5£,  a  function  of  x,  such  that  the  mean  squared  error  f 
is  a  minimum  where 

f  =  (y-x)T(y-x)  (2.8) 

and  x  is  that  minimizing  value  of  x.  The  optimal  estimate  x 
has  the  following  expression  [1,  p.  133]: 

x  =  Q  HT(HQHT+C)-1x  (2.9) 

M-l 

Let  fij=col  [tj  ,Xitj  ,  .  .  .  ,Xi  t_j  ]  ,  (l ,  j)  =1 ,  .  .  .  ,M 

t j=col [l,Xj,Xj,...,Aj  x] ,M=vector 

X  i=exp  [  /=TX  0  ( i - 1)  ],  X  Q= 2tt /M . 

—  —  2 

Let  and  c^ ^  ,  l<^i ,  j<M  denote  the  M  eigenvalues  of  the 

block  circulant  matrices,  H,  B(0)  and  C  respectively.  Then  from 

the  theory  of  circulant  matrices, 


where 


^.=001  [exp  (/-“l  ^j-(  (i-l)k+(  j-1  )0#  ( k  ,  v }  C  N  ] 

and  h=col(h.  . , 0< (i , j ) <p, (i , j ) ? (0 , 0) ] 

-  1  /  J 

and  p  is  the  width  of  the  PSF.  Using  the  fact  that  the  Fourier 
vectors  are  the  eigenvectors  of  Q,  H,  and  C,  with  the  corre- 

sponding  eigenvalues  p/|  j  y ^ j  j  |  ,  h^j  and  c^ ^ ,  the  following 

/N 

expression  can  be  written  for  x: 


x  = 


M 


n  h)^h 


h.  . 
11 


2+C  , 


ID1  1  ID 


2) ) f*T 
~1D 


(2.11) 


where  *  denotes  the  complex  conjugate  operator.  The  computation 

A 

of  x  can  be  done  by  using  FFT  and  the  computational  complexity 
2 

is  0(M  log  M) .  The  following  special  cases  of  (2.1)  are  of 
interest . 

Case  i  (Additive  noise  n  is  white  with  variance  y) : 

/s 

For  this  case,  C=y  I  and  x  reduces  to 


;  -  fljteK‘j/|p'lhi: 
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ID 
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T 
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(2.12) 


Case  ii  (Degradation  is  due  to  white  noise  only) : 
For  this  case,  H=I  and  hence  x  reduces  to 


x 


_1_  I  f  _ P 

M2  ~ij  (p+y|  ! Mi j 1  I2) 


x 


(2.13) 


The  restoration  algorithms  given  above  have  a  general  structure 
that  is  valid  for  the  so-called  causal,  semicausal  and  non-causal 
neighbor  sets.  (Actually,  due  to  the  torus  assumption,  the  neigh¬ 
bor  set  N={ (-1 , 0) , (0 , -1) , (-1 ,-l) }  does  not  correspond  to  a 


strictly  causal  model,  but  the  error  in  approximation  due  to  torus 
structure  is  very  small  for  causal  neighbor  sets.)  We  have  not 
made  any  special  assumptions  such  as  isotropy,  etc. 

Thus  far  we  have  assumed  that  the  parameters  0, p  characterizing 
the  original  image  are  known.  In  practice,  however,  they  are 
estimated  from  the  original  image  and  the  estimated  parameters 
are  used  in  place  of  true  parameters  in  the  restoration  algorithms. 


2 . 2  Parameter  Estimation 

The  popular  methods  of  estimation  are  the  method  of  least 
squares  (LS)  and  of  maximum  likelihood  (ML) .  However,  the 
classical  LS  estimates  are  not  in  general  consistent,  when  we 
consider  SAR  models  that  include  neighbor  set  dependence  in  all 
the  directions  [8] .  The  ML  estimation  scheme  involves  obtaining 
an  expression  for  the  likelihood  of  the  observations.  We  first 
assume  that  ( w ( - ) }  is  distributed  normally  with  zero  mean  and  unit 
variance.  To  ensure  stationarity ,  assume  |y^j(9)|^0,  i,j=l,...,M. 

By  using  the  property  of  the  noise  (w(*)}  and  (2.1),  the  log 
likelihood  of  observations  £np(yj9,p)  can  be  written  as 

«,n  p  (y  |  9  ,  P)  =  Hn  det  B (9)  -  (M2/2)  £n2irp  -  -L  ^ (y (s) -9Tz  (s) ) 2  (2.14) 

where 

z  (s)=col[y(s+(k,«.) )  ,  (k,Jt)tN] 

But  «.n  det  B(9)  =  ^(y^)  (2.15) 

Substitution  of  (2.15)  into  (2.14)  gives  the  desired  log  likeli¬ 
hood  function.  Due  to  the  log  likelihood  function  being  non¬ 
quadratic  in  9 ,  a  gradient  procedure  such  as  Newton-Raphson  may 
be  used.  Using  "off  the  shelf"  algorithms  may  be  computationally 
expensive.  We  use  the  iterative  scheme  given  below  in  (2.16)  and 
(2.17)  to  estimate  the  parameters: 


\+l=[R  -  S]_1(v  -  3-u) 

?t  ~  '  pt~ 

H  ■  jl  sL  <y<s)-et+l;(s>)2 

M 


(2.16) 

(2.17) 


w/  '  -■ 


?o  *  f’1"-  p0  -  ^  sen 

M 


v-  E  9ii'?=  E  (?ij?ij-?ii?ij) ' 

(i,j)£fi  13  (i,j)£fl  13  13  13  13 


C^j=col[[cos(^((i-l)k+(j-l)£))]  ,  (k,£)fcN]  ,  M-vector 


sT  .=col  [sin  (^-(  (i-1)  k+  ( j-l)£  ) )  ]  ,  (k,Jt)fcN]  ,  M-vector 
- 1]  M 


S  -  Ez(s)z(s)  ,  u=  E  z(s)y(s) 


Let  (0^,P^)  denote  the  final  estimates  of  (6,p).  In  the  implemen¬ 
tation  of  the  restoration  filters  the  estimates  9  and  p"  obtained 
from  Theorem  1  are  used  in  the  place  of  time  parameters.  For 
large  values  of  M,  the  estimates  S'  and  p-  are  close  to  the  ML 
estimates  and  are  asymptotically  consistent  and  efficient.  More 
discussion  regarding  the  estimation  scheme  may  be  found  in  [15-16] 


2 . 3  Error  Analysis  of  the  Restoration  Algorithms 

A  desirable  quantity  to  compute  is  the  expected  value  of  the 
minimum  error  .  Since  all  our  experiments  were  done  for  the  case 
of  degradation  due  to  additive  white  noise,  assume  from  now  on  that 
c  =  yl  ,  (2.18) 

where  y  is  the  noise  variance.  By  substitution  of  (2.9)  into 
(2.8),  the  following  expression  is  obtained  for  the  expected 
value  of  the  minimum  error  e  [1] : 

e  =  -4  Tr [Q-L  H  Q]  (2.19) 

M  ~  ~  ~  ~ 

where 

t  t  —  1 

L  =  Q  H  (H  Q  Ha+C) 


Using  the  fact  that  the  trace  of  a  matrix  is  the  sum  of  its 


eigenvalues, 
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e  =  — ^  Tr [Q-L  H  Q] 
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which  after  simplification  gives 
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Examples  of  Restoration 


The  algorithms  described  in  Section  2  were  implemented  on 
different  sizes  of  the  standard  USC  "girl  image"  under  different 
conditions . 

Example  1  (Image  with  additive  white  Gaussian  noise) 

Figure  1(a)  shows  the  original  girl  image  of  size  256x256 
and  intensity  variations  over  the  range  0-255.  Since  a  256x256 
image  may  not  be  stationary,  the  image  was  divided  into  16  blocks, 
each  of  size  64x64,  and  each  block  was  modeled  by  a  SAR  model, 
with  N={  (-1,0), (1,0), (0,-1), (0,1)}.  Parameters  were  estimated 
using  the  iterative  scheme  in  Section  2.2.  Figure  1(b)  shows 
the  noisy  image  with  SNR=7  db.  The  restoration  algorithm  in  (2.13) 
was  implemented  and  each  block  was  restored  independently.  We 
emphasize  that  blocking  was  done  due  to  considerations  of 
stationarity  alone  and  not  due  to  any  computational  problems. 

The  restored  image  is  shown  in  Figure  1(c). 

Due  to  the  blocking,  artificial  lines  are  present  in  the 
restored  image  at  known  locations.  These  lines  may  be  removed 
either  in  the  Fourier  transform  domain  or  by  a  linear  inter¬ 
polation  technique  (2J.  We  used  the  latter  technique.  To  remove 
vertical  lines  (3  pixels  in  width),  each  pixel  in  the  lines  was 
replaced  by  a  linear  weighted  sum  of  neighboring  pixels  (including 
the  current  pixel)  with  weights  0.25,  0.2,  0.1,  0.2,  0.25.  The 
image  with  vertical  bars  removed  is  in  Figure  1(d).  The  same 
technique  was  used  to  remove  the  horizontal  bars  and  the  final 
restored  picture  is  shown  in  Figure  1 (e) . 
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Example  2  (Image  with  blur  and  noise) : 

For  this  experiment,  a  64x64  window  of  the  original  girl's 
face  was  used.  The  parameters  corresponding  to  different  neighbor 
sets  were  computed.  The  image  was  blurred  by  using  the  PSF 

h(k,M  =  ^  exp{-0.4(k2+£2) }  (3.1) 

and  Gaussian  noise  of  SNR=7  db  was  added.  The  results  of 
implementing  (2.12)  are  given  in  Figure  2  for  different  neighbor 
sets  of  dependence.  Similar  results  for  SNR=0  db  are  shown  in 
Figure  3.  The  details  of  the  neighbor  sets  used  may  be  found  in 
the  figure  caption. 

To  get  some  idea  as  to  the  amount  of  theoretical  improvement 

that  is  possible,  when  using  different  neighbor  sets,  numerical  values 

of  the  expected  error  e  were  computed  using  (2.20)  for  those  neighbor 

sets.  The  quantity  g  defined  as 

g  =  10  log  MSE  between  original  and  degraded 

Q 

is  tabulated  in  Table  1  for  SNR=7  db  and  0  db  and  for  the  blur 
function  with  PSF  given  in  (3.1).  The  values  of  the  numerator  for 
SNR=7  db  and  0  db  were  511.41  and  1910.1  respectively.  Although 
not  reported  here,  the  implementation  for  the  case  of  colored 
noise  can  be  done  with  equal  ease. 


Table  1 


Performance  bounds  in  decibels  as  predicted  by  various  SAR 
models  for  the  PSF  in  (3.1)  and  M=64. 


Neighbor  set  N  SNR=7  db  SNR: 

{  (-1,0)  ,  (0,-1)  ,  (-1,-1)  }  causal  4.687  8 

{ (-1,0) , (1,0) , (0,-1) Jsemi-causal  6.020  9 

{ (-1 , 0) , ( 1 , 0) , ( 0 , -1) , (0 , 1) }  non-causal  5.634  8 

( (-1,0) ,  (0,1)  ,  (1,1)  , (1,0) , (0,-1) , (-1,-1)  } 

non-causal  5.530  8 

{  (-1,0)  , (-1,1)  , (0,1)  , (1,0)  ,  (1,-1) , (0,-1) } 

non-causal  4.630  6 

{  (-1,0)  ,  (-1,1) , (0,1) , (1,1) , (1,0) , (1,-1)  , 

(0 , -1) , (-1 , -1) }  non-causal  5.272  12 


0  db 
678 
955 
641 

531 

811 

12 


4.  Discussion 


We  have  developed  restoration  schemes  using  spatial  auto¬ 
regressive  random  field  models.  When  dependence  in  all  directions 
is  included,  the  neighbors  corresponding  to  the  boundary  pixels 
are  not  defined  for  these  models.  Hence  some  assumptions  must  be 
made  regarding  the  distribution  of  these  boundary  pixels.  It  would 
be  beneficial  to  make  assumptions  which  aid  in  reducing  excessive 
computations  otherwise  present  in  MMSE  restoration  problems,  while 
at  the  same  time  not  sacrificing  any  optimality  properties.  This 
leads  to  the  problem  of  finding  covariance  matrices  of  SAR  models 
that  are  diagonalized  by  fast  transforms  such  as  the  FFT ,  sine  or 
cosine  transforms.  By  appropriately  imposing  some  assumptions 
about  the  boundary  pixels,  covariance  matrices  which  are  diagonalized 
by  these  fast  trasnforms  can  be  realized  [10] .  We  have  used  one 
such  assumption,  viz.,  that  the  images  are  represented  on  a 
toroidal  lattice.  The  use  of  the  torus  lattice  representation 
leads  to  the  use  of  the  FFT  in  the  implementation  of  the  filters, 
for  any  arbitrary  neighbor  set. 

Such  fast  restoration  filters  have  been  considered  in  the 
literature  [6-7]  for  different  neighbor  sets,  using  fast  sine 
transforms.  There  are  several  disadvantages  in  using  this  formu¬ 
lation.  Since  the  fast  sine  trasnforms  diagonalize  symmetric, 
tridiagonal  Toeplitz  matrices,  the  decomposition  schemes  suggested 
in  [7]  are  valid  only  for  isotropic  neighbor  sets.  For  instance, 
for  the  neighbor  set  N  =  {  (0  ,-l)  ,  (1 ,0)  ,  (-1 , 0) },  it  is  required 
that  0.  n=0_1  _.  However,  our  experiments  with  this  model  and 


the  estimation  scheme  in  Section  2  do  not  support  this  assumption. 
Also,  when  neighbor  sets  like  {  (-1 , 0) , (1 , 0) , (0 , -1) ,  (0 , 1)  ,  (1 , 1)  , 
(-1,-1) , (-1,1) , (1,-1) , (0,2) , (0,-2) , (-2,0) , (2,0) }  are  considered, 
the  resulting  covariance  matrices  have  banded  Toeplitz  structures 
and  hence  are  not  diagonalized  exactly  by  sine  transforms,  so  that 
approximate  methods  are  required.  On  the  other  hand,  the  approach 
taken  in  this  paper  does  not  involve  any  such  approximations.  Given 
that  some  assumptions  have  to  be  made  regarding  the  boundary 
pixels,  it  is  preferable  to  make  assumptions  which  have  more 
general  applicability. 

The  assumption  of  torus  lattices  can  also  be  justified  as 
follows:  We  computed  the  noramlized  autocorrelation  function  at 
lower  lags  for  different  SAR  models  using  the  torus  assumption 
and  also  by  Fourier  inversion  of  the  spectral  density  function 
corresponding  to  some  SAR  models  without  the  torus  assumption. 

The  numerical  values  differ  in  the  second  or  third  decimal  place. 
Also,  the  torus  assumption  leads  to  an  intuitively  appealing 
relationship,  viz.,  the  eigenvalues  of  the  covariance  matrix  are 
the  two-dimensional  discrete  spectral  density  function. 

We  have  also  estimated  the  parameters  of  the  underlying  SAR 
models.  The  need  for  good  estimates  can  be  explained  as  follows. 
There  is  always  a  finite  error  in  restoration;  this  error  will  be 
greater  if  the  parameters  characterizing  the  SAR  models  are  not 
consistent  or  efficient.  Though  estimation  schemes  for  unilateral 
neighbor  sets  are  simple,  they  are  not  so  for  the  SAR  models 
considered  here.  By  using  the  asymptotically  consistent  and 


r 


efficient  estimates  considered  in  this  paper,  the  errors  due  to 

incorrect  model  specification  can  bo  reduced. 

Fast  Fourier  transforms  for  image  restoration  using  Wiener 

filters  [1,17]  have  been  considered  earlier  by  approximating  the 

block  Toeplitz  covariance  structures  by  block  circulant  matrices. 

The  important  difference  between  the  scheme  developed  here  and 

the  earlier  schemes  is  in  the  characterization  of  signal  statistics. 

The  MMSE  schemes  need  exact  knowledge  of  the  two-dimensional 

spectral  density  function  (SDF) .  In  practice  it  is  not  known. 

Hence,  the  true  SDF  is  replaced  by  an  estimate.  Typical  estimates 

of  the  SDF  can  be  obtained  by  making  suitable  assumptions  regarding 

the  ACF  as  in  [17] ,  or  by  using  estimates  such  as  the  FFTs  of 

empirical  covarance  matrices.  The  former  scheme  assumes  very 

simple  functions  for  the  ACF,  while  the  latter  uses  inconsistent 

estimates  of  the  SDF.  In  our  approach,  the  true  two-dimensional 

SDF  is  an  explicit  function  of  the  parameters  given  by 

2 

sy  (?'p,vl'v2)  =  ,  where  v1=XQi,  v2=XQ  j .  By  using 

asymptotically  consistent  estimates  like  S^  (F,"p, , v^)  ,  as  done 
here,  better  performance  can  be  achieved.  Also,  the  underlying 
structures  of  SAR  models  are  more  complicated  and  varied  than  that 
of  the  simple  causal  model  in  [17] . 
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b)  noisy  SNR  =  7  db 

c)  restored  image  (note  the  presence 
of  lines  due  to  blocking) 

d)  vertical  lines  removed 


e)  final  restored  image. 
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Restoration  of  image  '  '  _ 

containing  blur  (4.1) 
and  Gaussian  noise 
Original ; 

Blurred  image; 

Restored  image,  SNR  =  a; 

Blur  and  noise  (7  db) ;  , 

Restoration  with  N  =  {  (0,-1)  ,  (  1,0),  (  1,-1)  ), 

N  =  {  (-1,0)  (1,0)  ,  (0,-D  }; 

N  =  {  (0,1)  , (0,-1)  , (-1,0), (1,0) } ; 

N  =  M0, 1), (1,0), (0,-D, (-1,0), (0,-1  ,  H,  I- 
N  =  {  (0,1)  ,  (0,-1)  ,  (1,0)  ,  (1,1)  ,  (-1,-1)  ,  (-1-,1)  ,  (1/  Iw* 


Similar  to  Fig.  2  with 
SNR  =  0  db. 


